#ifndef _MAGTET_H
#define _MAGTET_H

#include "gctl/core.h"
#include "gctl/geometry.h"
#include "gctl/io.h"
#include "gctl/utility.h"

#include "omp.h"

#define Cm 1e-7
#define TOL 1e-10

#define TETFILE  "tet-file"
#define MAGFILE  "mag-file"
#define SITEFILE "site-file"
#define OBSFILE  "obs-file"
#define CALTYPE  "cal-type" // potential, gradient or tensor

#define RUN_ECHO(action, msg) do {std::clog << msg << "... \n"; action;} while(0);

struct magtet_para
{
	double mag_amp[4];
	gctl::point3dc fnorm[4], enorm[12], etang[12];
};

typedef gctl::type_tetrahedron<magtet_para> mag_tetrahedron;

class magtet
{
public:
	magtet(){}
	virtual ~magtet(){}

	void read_tet(std::string tet_name);
	void init_magz(std::string para);
	void init_site(std::string para);
	void write_text(std::string out_name);

	void cal_tensors();
	double mag_potential(const mag_tetrahedron &tet, const gctl::point3dc &site);
	gctl::point3dc mag_gradient(const mag_tetrahedron &tet, const gctl::point3dc &site);
	gctl::tensor mag_tensor(const mag_tetrahedron &tet, const gctl::point3dc &site);
	void total_potential();
	void total_gradient();
	void total_tensor();

	//double sign(double s);
	void routine(const char *para_file);

protected:
	int node_num_, ele_num_, site_num_;
	gctl::array<gctl::vertex3dc> node_; // 四面体元素的顶点集
	gctl::array<mag_tetrahedron> ele_; // 四面体的元素集
	gctl::array<gctl::point3dc> magz_; // 磁化矢量数组
	gctl::array<magtet_para> ele_para_; // 四面体的张量属性与磁化矢量
	gctl::array<gctl::point3dc> site_; // 观测点位置集

	gctl::array<double> mag_pot_; // 正演磁位
	gctl::array<gctl::point3dc> mag_grad_; // 正演磁场分量
	gctl::array<gctl::tensor> mag_tensor_; // 正演磁场张量

	gctl::getoption gopt_; // 参数文件解释对象
};

#endif // _MAGTET_H